# Set working directory to main replication folder
setwd("C:/Users/Eric/Dropbox/Book/analysis/replication/")

# Load packages
library(data.table)
library(msm)
library(plyr)
library(dplyr)
library(tidyverse)
library(survival)
library(ggplot2)
library(stargazer)
library(strucchange)
library(xtable)

# Load functions
source("./function_code/msm_coef_plot.R")
source("./function_code/msm_coef_table.R")
source("./function_code/hazard_coef_plot.R")

# Read in main dataset
d = fread("./data/wardata_book.csv")
d = d |> arrange(warName, dayNum)

# Read in data merged with ICM dataset
d3 = fread("./data/tpidata_book.csv")

# Negotiation start data
s = fread("./data/negstarts.csv")


#### THE VALIDITY OF THE 1945 LINE ####

## Structural break test ##
# Create base data
yr = 1820:2005 
nw = nwd = nnd = NA
for (i in 1:length(yr)) {
  oy = d[format(d$date, "%Y")==yr[i],]
  nw[i] = length(unique(oy$warName))   # Number of wars
  nwd[i] = nrow(oy)                    # Number of war-days
  nnd[i] = sum(oy$neg)                 # Number of negotiation days
}

# Create time series
negs_ts = ts(nnd, start=1820)

# Run the break point model
breaktest = breakpoints(negs_ts ~ nw + nwd)

# Identify break points
#### TABLE A3: Results of structural break point models ####
summary(breaktest) # Best model in BIC has two cutpoints at 1945 and 1972

rm(oy, yr, nw, nwd, nnd, negs_ts)

## Logistic regressions ##
y = 1824:2002
bics = est = se = NA
for (i in 1:length(y)) {
  # Run logistic regression using year Y as the threshold for a binary variable
  ntest = glm(neg ~ I(format(date, "%Y") > y[i]), data=d, family=binomial)
  bics[i] = BIC(ntest)
  est[i] = summary(ntest)$coefficients[2,1]
  se[i] = summary(ntest)$coefficients[2,2]
}
bicout = data.frame(year=y, bic=bics, est, se)

# Identify year(s) with lowest BIC
bicout$year[which(bicout$bic==min(bicout$bic))] # 1945 and 1946

rm(y, bics, est, se)


#### MODELING WARTIME NEGOTIATIONS ####

## Multistate models ##
d$negState3 = factor(d$negState3)

# State transition table
Qs = statetable.msm(negState3, warName, data=d)

# Make the crude initial values
Qinit3 = matrix(0, ncol=3, nrow=3)
for (i in 1:nrow(Qs)) {
  rowtotal = sum(Qs[i,])
  for (j in 1:ncol(Qs)) {
    Qinit3[i,j] = Qs[i,j]/rowtotal 
  }
}

# Multistate model for pre-1945 wars
threeFitPre = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==0), qmatrix=Qinit3,
                  covariates = ~ AOS60 + AOA + issues + contig + 
                    cincP + polBinA + keyAllies + nBellig + logBattles + anyDipEnemy, 
                  control=list(fnscale=4000, maxit=500))

# Multistate model for post-1945 wars
threeFitPost = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==1), qmatrix=Qinit3,
                   covariates = ~ AOS60 + AOA + issues + contig + 
                     cincP + polBinA + nuk + keyAllies + nBellig + logBattles + anyDipEnemy, 
                   control=list(fnscale=4000, maxit=2000))

#### FIGURE 3.2: Coefficient plots for multistate models of negotiation in pre-1945 wars (n = 23,711) ####
threeFitPre_res = msmCoef(threeFitPre)
msmCoefPlot(threeFitPre_res, -4, 2.8, "Pre60") 
file.rename("./figures/msmCoefPlotPre60.pdf",
            "./figures/figure_3.2.pdf")

#### FIGURE 3.3: Coefficient plots for multistate models of negotiation in post-1945 wars (n = 13,123) ####
threeFitPost_res = msmCoef(threeFitPost)
msmCoefPlot(threeFitPost_res, -8, 6, "Post60") 
file.rename("./figures/msmCoefPlotPost60.pdf",
            "./figures/figure_3.3.pdf")


#### TABLE 3.1: Hazard ratios ####

# Pre-1945, start negotiations
round(hazard.msm(threeFitPre)[[1]][1,c(2,1,3)], 3)  # 1.068 1.151 1.241
  
# Pre-1945, negotiated settlement
round(hazard.msm(threeFitPre)[[1]][4,c(2,1,3)], 3)  # 1.114 1.269 1.447
  
# Post-1945, start negotiations
round(hazard.msm(threeFitPost)[[1]][1,c(2,1,3)], 3) # 0.883 1.092 1.350 

# Post-1945, negotiated settlement
round(hazard.msm(threeFitPost)[[1]][4,c(2,1,3)], 3) # 1.181 2.216 4.158 



#### Issue salience across eras ####
iss = d |> group_by(warName) |> dplyr::select(issues, post45) |> slice(1)
round(mean(iss$issues[iss$post45==0]), 3) # 1.877 
round(mean(iss$issues[iss$post45==1]), 3) # 2.229
t.test(iss$issues[iss$post45==0], iss$issues[iss$post45==1])  # p = 0.08073
rm(iss)


#### FIGURE 3.4: Smoothed plots of propensity to negotiate over the course of wars ####
wn = unique(d$warName)
rseq = seq(0, 1, 0.05)

warsPre = unique(d$warName[which(d$post45==0)])
warsPost = unique(d$warName[which(d$post45==1)])

# Calculate pre-1945 war trajectories
meansPre = list()
for (i in 1:length(warsPre)) {
  oneWar = d[d$warName==warsPre[i],]
  
  mneg = mnegInt = mnegExt = NA
  for (j in 1:length(rseq)) {
    part = oneWar[oneWar$propDoneR==rseq[j],]
    mneg[j] = mean(part$neg, na.rm=T)
    mnegInt[j] = mean(part$negInt, na.rm=T)
    mnegExt[j] = mean(part$negExt, na.rm=T)
  }
  meansPre[[i]] = data.frame(propDone=rseq, 
                             neg=mneg, negInt=mnegInt, negExt=mnegExt)
}
mrPre = do.call('rbind', meansPre)
mrPre$neg = na.approx(mrPre$neg, na.rm=F)
mrPre$negInt = na.approx(mrPre$negInt, na.rm=F)
mrPre$negExt = na.approx(mrPre$negExt, na.rm=F)
mrPre = na.omit(mrPre)
smPre = ksmooth(mrPre$propDone, mrPre$neg, bandwidth=0.2)
smPreInt = ksmooth(mrPre$propDone, mrPre$negInt, bandwidth=0.2)
smPreExt = ksmooth(mrPre$propDone, mrPre$negExt, bandwidth=0.2)

trajPre = data.frame(propDone=smPre$x, neg=smPre$y, negInt=smPreInt$y, negExt=smPreExt$y)
trajPre$propDone = round_any(trajPre$propDone, 0.05, floor)
trajPre = unique(trajPre)
trajPre = trajPre[!duplicated(trajPre$propDone),]
trajPre$era = "Pre-1945"

# Calculate post-1945 war trajectories
meansPost = list()
for (i in 1:length(warsPost)) {
  oneWar = d[d$warName==warsPost[i],]
  
  mneg = mnegInt = mnegExt = NA
  for (j in 1:length(rseq)) {
    part = oneWar[oneWar$propDoneR==rseq[j],]
    mneg[j] = mean(part$neg, na.rm=T)
    mnegInt[j] = mean(part$negInt, na.rm=T)
    mnegExt[j] = mean(part$negExt, na.rm=T)
  }
  meansPost[[i]] = data.frame(propDone=rseq, 
                              neg=mneg, negInt=mnegInt, negExt=mnegExt)
}
mrPost = do.call('rbind', meansPost)
mrPost$neg = na.approx(mrPost$neg, na.rm=F)
mrPost$negInt = na.approx(mrPost$negInt, na.rm=F)
mrPost$negExt = na.approx(mrPost$negExt, na.rm=F)
mrPost = na.omit(mrPost)
smPost = ksmooth(mrPost$propDone, mrPost$neg, bandwidth=0.2)
smPostInt = ksmooth(mrPost$propDone, mrPost$negInt, bandwidth=0.2)
smPostExt = ksmooth(mrPost$propDone, mrPost$negExt, bandwidth=0.2)

trajPost = data.frame(propDone=smPost$x, neg=smPost$y, negInt=smPostInt$y, negExt=smPostExt$y)
trajPost$propDone = round_any(trajPost$propDone, 0.05, floor)
trajPost = unique(trajPost)
trajPost = trajPost[!duplicated(trajPost$propDone),]
trajPost$era = "Post-1945"

# Put pre-1945 and post-1945 trajectory datasets together
traj = rbind(trajPre, trajPost); rm(trajPre, trajPost)
traj$era = factor(traj$era, levels=c("Pre-1945", "Post-1945"))

trajPlot = traj |> pivot_longer(cols=contains("neg"), names_to="typeNeg", values_to="propNeg") 
trajPlot$typeNeg = factor(trajPlot$typeNeg, levels=c("neg", "negInt", "negExt"))

# Make and export the plot
ggplot(trajPlot, aes(x=propDone, y=propNeg)) + geom_line(aes(color=typeNeg), lineend="round", linewidth=1.5) + theme_bw() + 
  xlab("Proportion of war elapsed") + ylab("Pr(Negotiation)") +
  facet_grid(rows=vars(era)) +
  scale_color_manual("Negotiation", values=c("black", "gray35", "gray70"), 
                     labels=c("All", "Internal", "External")) + 
  theme(legend.position="bottom",
        legend.direction = "horizontal") +
  ylim(0,0.45) 
ggsave("./figures/figure_3.4.pdf", height=6.5, width=4.5)


#### MANIFESTATIONS OF PRESSURE ####

# Set up data
survWar = Surv(as.numeric(d$dayNum), as.numeric(d$dayNumNext), as.numeric(d$terminate))
controls1 = "AOA + issues + contig + cincP + polBinA + nuk + nBellig + anyDipEnemy + keyAllies + logBattles"

# Model with all negotiations
endFit0a = coxph(as.formula(paste("survWar ~ neg + AOS60 +  
                                  logBattles:logDayNum + AOA:logDayNum", controls1, sep="+")), cluster=warName, data=d)
summary(endFit0a)
cox.zph(endFit0a)

# Impact of negotiation on war termination 
round(coef(endFit0a)['neg'], 3)         # Raw: 1.977
round(exp(coef(endFit0a)['neg']), 3)    # Hazard: 7.221


# Model that disaggregates internal and external negotiations
endFit0b = coxph(as.formula(paste("survWar ~ negInt + negExt + AOS60 +  
                                  logBattles:logDayNum + AOA:logDayNum", controls1, sep="+")), cluster=warName, data=d)
summary(endFit0b)
cox.zph(endFit0b)

# Impact of negotiation on war termination
round(coef(endFit0b)['negInt'], 3)      # Internal, raw: 2.466
round(exp(coef(endFit0b)['negInt']), 3) # Internal, hazard: 11.779

round(coef(endFit0b)['negExt'], 3)      # External, raw: 1.045
round(exp(coef(endFit0b)['negExt']), 3) # External, hazard: 2.843

# Linear hypothesis test for whether difference in 
# internal and external negotiations are statistically significant
car::linearHypothesis(endFit0b, "negInt=negExt") # p = 0.00076


#### FIGURE 3.5: Hazard models of war termination (n = 36,834) ####
hazardCoefPlot(model1=endFit0a, model2=endFit0b, fileID="Main")
file.rename("./figures/hazardCoefPlotMain.pdf",
            "./figures/figure_3.5.pdf")


#### MINOR SOURCES OF PRESSURE ####

# Minor power designations by war
minors = d |> dplyr::group_by(warName) |> dplyr::summarize(minor=median(minor), post45=median(post45))
p45_minor = table(minors$post45, minors$minor)
chisq.test(p45_minor) # p = 0.744

# Are pressured negotiations associated with minor powers?
press_minor = merge(s, minors, by="warName")
table(press_minor$minor, press_minor$pressure)
round(prop.table(table(press_minor$minor, press_minor$pressure), 1), 3)


# Make initialization matrix
Qminor = statetable.msm(negStateIE, warName, data=d)
QinitIE = matrix(0, ncol=3, nrow=3)
for (i in 1:nrow(Qminor)) {
  rowtotal = sum(Qminor[i,])
  for (j in 1:ncol(Qminor)) {
    QinitIE[i,j] = Qminor[i,j]/rowtotal 
  }
}

minorFit1a = msm(negStateIE ~ dayNum, subject=warName, data=d, qmatrix=QinitIE,   # ADD BACK 1945?
                 covariates = ~ minor + issues + contig + cincP + polBinA + nuk +
                   nBellig + keyAllies + logBattles + anyDipEnemy + AOS60 + AOA, 
                 control=list(fnscale=4000, maxit=500))
hazard.msm(minorFit1a)


#### FIGURE 3.6: Coefficient plots for multistate models of negotiation onset, disaggregated by third-party pressure (n = 36,864) ####
minorFit1a_res = msmCoef(minorFit1a)
msmCoefPlot2(minorFit1a_res, -1.5, 2.8, "SplitMinor") 
file.rename("./figures/msmCoefPlotSplitMinor.pdf",
            "./figures/figure_3.6.pdf")

# The impact of minor power wars on diplomatic intervention attempts
round(hazard.msm(minorFit1a)[["minor"]][1,1], 3)       # 0.709
round(exp(hazard.msm(minorFit1a)[["minor"]][1,1]), 3)  # 2.031

# Total interventions in the dataset
length(which(d3$intervention==1)) # 151

# Almost half (75) of all interventions occur in the first third of conflicts' durations
length(which(d3$intervention==1 & d3$propDone < 0.333))

# 61 of the 75 diplomatic efforts that take place in this time fail to get anywhere
table(d3$tpi_outcome_2[which(d3$intervention==1 & d3$propDone < 0.333)])


## What factors account for when third parties apply more pressure in post-1945 wars? ##
surv3dip = Surv(as.numeric(d3$dayNum), as.numeric(d3$dayNumNext), as.numeric(d3$intervention))

tpiFit0a = coxph(surv3dip ~ AOS60 + AOA + issues + contig + cincP + polBinA + nuk + minor + 
                   nBellig + keyAllies + logBattles + AOS60:logDayNum + strata(anyDipEnemy), cluster=warName, data=d3)
summary(tpiFit0a)
cox.zph(tpiFit0a)


# The impact of minor power wars on diplomatic intervention attempts
round(coef(tpiFit0a)['minor'], 3)      # 1.999
round(exp(coef(tpiFit0a)['minor']), 3) # 7.385

#### FIGURE 3.7: Coefficient plots for survival model of third-party diplomatic interventions in post-1945 wars (n = 11,398) ####
hazardCoefPlot(model1=tpiFit0a, model2="none", fileID="TPI")
file.rename("./figures/hazardCoefPlotTPI.pdf",
            "./figures/figure_3.7.pdf")

# What is the correlation between overall imbalance and the length of war?
cor(d$AOA, d$dayNum) # 0.531
cor(d$AOA[which(d$post45==1)], d$dayNum[which(d$post45==1)]) # 0.770


######################################################################################################


#### APPENDIX MATERIALS FOR CHAPTER 3 ####

##### APPENDIX 3.1: THE 1945 LINE #####

###### TABLE A3: Results of structural break point models ######

# See "The Validity of the 1945 Line" section above.

###### FIGURE A2: BIC measures for a series of logistic regressions splitting the negotiation data by year ######

bicout |> ggplot(aes(year, bic)) + geom_line() + theme_bw() +
  xlab("Year") + ylab("BIC")
ggsave("./figures/figure_a2.pdf", width=6, height=3.5)


##### APPENDIX 3.2: MULTISTATE MODELS #####

###### APPENDIX 3.2.1: Full 60-Day Results ######

####### TABLE A4: Results of multistate models for pre-1945 wars, using 60-day recent imbalance measures #######
msmCoefTable(threeFitPre, checkMinor = 0)


####### TABLE A5: Results of multistate models for post-1945 wars, using 60-day recent imbalance measures #######
msmCoefTable(threeFitPost, checkMinor = 0)


###### APPENDIX 3.2.2: Full 30-Day Results ######

# Pre-1945 wars
threeFitPre30 = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==0), qmatrix=Qinit3,
                  covariates = ~ AOS30 + AOA + issues + contig + 
                    cincP + polBinA + keyAllies + nBellig + logBattles + anyDipEnemy, 
                  control=list(fnscale=4000, maxit=500))
# Post-1945 wars
threeFitPost30 = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==1), qmatrix=Qinit3,
                   covariates = ~ AOS30 + AOA + issues + contig + 
                     cincP + polBinA + nuk + keyAllies + nBellig + logBattles + anyDipEnemy, 
                   control=list(fnscale=4000, maxit=2000))

####### TABLE A6: Results of multistate models for pre-1945 wars, using 30-day recent imbalance measure #######
threeFitPre30_res = msmCoef(threeFitPre30)
msmCoefTable(threeFitPre30, checkMinor = 0)

####### FIGURE A3: Coefficient plots for multistate models of negotiation in pre-1945 wars, using 30-day recent imbalance measure #######
msmCoefPlot(threeFitPre30_res, -4, 3, "Pre30") 
file.rename("./figures/msmCoefPlotPre30.pdf",
            "./figures/figure_a3.pdf")

####### TABLE A7: Results of multistate models for post-1945 wars, using 30-day recent imbalance measure #######
threeFitPost30_res = msmCoef(threeFitPost30)
msmCoefTable(threeFitPost30, checkMinor = 0)

####### FIGURE A4: Coefficient plots for multistate models of negotiation in post-1945 wars, using 30-day recent imbalance measure #######
msmCoefPlot(threeFitPost30_res, -8, 6, "Post30") 
file.rename("./figures/msmCoefPlotPost30.pdf",
            "./figures/figure_a4.pdf")


###### APPENDIX 3.2.3: Full 90-Day Results ######

# Pre-1945 wars
threeFitPre90 = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==0), qmatrix=Qinit3,
                    covariates = ~ AOS90 + AOA + issues + contig + 
                      cincP + polBinA + keyAllies + nBellig + logBattles + anyDipEnemy, 
                    control=list(fnscale=4000, maxit=500))
# Post-1945 wars
threeFitPost90 = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==1), qmatrix=Qinit3,
                     covariates = ~ AOS90 + AOA + issues + contig + 
                       cincP + polBinA + nuk + keyAllies + nBellig + logBattles + anyDipEnemy, 
                     control=list(fnscale=4000, maxit=2000))

####### TABLE A8: Results of multistate models for pre-1945 wars, using 90-day recent imbalance measure #######
threeFitPre90_res = msmCoef(threeFitPre90)
msmCoefTable(threeFitPre90, checkMinor = 0)

####### FIGURE A5: Coefficient plots for multistate models of negotiation in pre-1945 wars, using 90-day recent imbalance measure #######
msmCoefPlot(threeFitPre90_res, -4, 3, "Pre90") 
file.rename("./figures/msmCoefPlotPre90.pdf",
            "./figures/figure_a5.pdf")

####### TABLE A9: Results of multistate models for pre-1945 wars, using 90-day recent imbalance measure #######
threeFitPost90_res = msmCoef(threeFitPost90)
msmCoefTable(threeFitPost90, checkMinor = 0)

####### FIGURE A6: Coefficient plots for multistate models of negotiation in post-1945 wars, using 90-day recent imbalance measure #######
msmCoefPlot(threeFitPost90_res, -8, 5, "Post90") 
file.rename("./figures/msmCoefPlotPost90.pdf",
            "./figures/figure_a6.pdf")

#### APPENDIX 3.3: HAZARD MODELS ####

##### TABLE A10: Results of Cox proportional hazard models on war termination #####
stargazer(endFit0a, endFit0b,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum"),
          order = c("neg", "negInt", "negExt",
                    "issues", "contig", "cincP", "polBinA", "nuk",
                    "nBellig", "anyDipEnemy", "keyAllies", "logBattles", 
                    "AOS60", "AOA"),
          covariate.labels = c("Negotiation", "Internal negotiation", "External negotiation", 
                               "Issue salience", "Contiguity",
                               "CINC ratio", "Democratic initiator", "Nuclear",
                               "Number of states", 
                               "Opp. diplomatic representation", 
                               "Major allies",
                               "Completed battles", 
                               "Recent imbalance", "Overall imbalance"),
          dep.var.labels.include=T,
          dep.var.labels = c("War termination", "War termination")
)

#### APPENDIX 3.4: ALTERNATIVE MEASURE OF CAPABILITIES ####

# Pre-1945 wars
threeFitPreDOE = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==0), qmatrix=Qinit3,
                    covariates = ~ AOS60 + AOA + issues + contig + 
                      doeWinLose + 
                      polBinA + nBellig + keyAllies + logBattles + anyDipEnemy, 
                    control=list(fnscale=4000, maxit=500))

# Post-1945 wars
threeFitPostDOE = msm(negState3 ~ dayNum, subject=warName, data=subset(d, post45==1), qmatrix=Qinit3,
                     covariates = ~ AOS60 + AOA + issues + contig + 
                       doeWinLose + 
                       polBinA + nuk + nBellig + keyAllies + logBattles + anyDipEnemy, 
                     control=list(fnscale=5000, maxit=5000))

#### TABLE A11: Results of multistate models for pre-1945 wars, replacing CINC ratios with DOE scores ####
threeFitPreDOE_res = msmCoef(threeFitPreDOE)
msmCoefTable(threeFitPreDOE, checkMinor = 0)

#### FIGURE A7: Coefficient plots for multistate models of negotiation in pre-1945 wars, replacing CINC ratios with DOE scores ####
msmCoefPlot(threeFitPreDOE_res, -3, 4, "PreDOE") 
file.rename("./figures/msmCoefPlotPreDOE.pdf",
            "./figures/figure_a7.pdf")

#### TABLE A12: Results of multistate models for post-1945 wars, replacing CINC ratios with DOE scores ####
threeFitPostDOE_res = msmCoef(threeFitPostDOE)
msmCoefTable(threeFitPostDOE, checkMinor = 0)

#### FIGURE A8: Coefficient plots for multistate models of negotiation in post-1945 wars, replacing CINC ratios with DOE scores ####
msmCoefPlot(threeFitPostDOE_res, -65, 70, "PostDOE") 
file.rename("./figures/msmCoefPlotPostDOE.pdf",
            "./figures/figure_a8.pdf")

## DOE results for hazard model of war termination ##
controls2 = "AOA + issues + contig + doeWinLose + polBinA + nuk + nBellig + keyAllies + logBattles + anyDipEnemy"

# Model with all negotiations
endFitDOEa = coxph(as.formula(paste("survWar ~ neg + AOS60 +  
                                  logBattles:logDayNum + AOA:logDayNum", controls2, sep="+")), cluster=warName, data=d)
cox.zph(endFitDOEa)

# Impact of negotiation on war termination 
round(coef(endFitDOEa)['neg'], 3)         # Raw: 1.958
round(exp(coef(endFitDOEa)['neg']), 3)    # Hazard: 7.084


# Model that disaggregates internal and external negotiations
endFitDOEb = coxph(as.formula(paste("survWar ~ negInt + negExt + AOS60 +  
                                  logBattles:logDayNum + AOA:logDayNum", controls2, sep="+")), cluster=warName, data=d)
cox.zph(endFitDOEb)

# Impact of negotiation on war termination
round(coef(endFitDOEb)['negInt'], 3)      # Internal, raw: 2.463
round(exp(coef(endFitDOEb)['negInt']), 3) # Internal, hazard: 11.739 

round(coef(endFitDOEb)['negExt'], 3)      # External, raw: 1.031
round(exp(coef(endFitDOEb)['negExt']), 3) # External, hazard: 2.803 

# Linear hypothesis test for whether difference in 
# internal and external negotiations are statistically significant
car::linearHypothesis(endFitDOEb, "negInt=negExt") # p = 0.0005


#### TABLE A13: Results of Cox proportional hazard models on war termination, replacing CINC ratios with DOE scores ####
stargazer(endFitDOEa, endFitDOEb,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum"),
          order = c("neg", "negInt", "negExt",
                    "issues", "contig", "doeWinLose", "polBinA", "nuk",
                    "nBellig", "anyDipEnemy", "keyAllies", "logBattles", 
                    "AOS60", "AOA"),
          covariate.labels = c("Negotiation", "Internal negotiation", "External negotiation", 
                               "Issue salience", "Contiguity",
                               "DOE score", "Democratic initiator", "Nuclear",
                               "Number of states", 
                               "Opp. diplomatic representation", 
                               "Major allies",
                               "Completed battles", 
                               "Recent imbalance", "Overall imbalance"),
          dep.var.labels.include=T,
          dep.var.labels = c("War termination", "War termination")
)

#### FIGURE A9: Hazard models of war termination, replacing CINC ratios with DOE scores ####
hazardCoefPlot(model1=endFitDOEa, model2=endFitDOEb, fileID="DOE")
file.rename("./figures/hazardCoefPlotDOE.pdf",
            "./figures/figure_a9.pdf")


#### APPENDIX 3.5: MINOR POWER WARS ####

##### TABLE A14: Results of multistate models on minor power wars #####
msmCoefTable(minorFit1a, checkMinor = 1)


#### APPENDIX 3.6: THIRD-PARTY DIPLOMATIC INTERVENTION ####

##### TABLE A15: Results of a Cox proportional hazard model on third-party diplomatic intervention in post-1945 wars #####
stargazer(tpiFit0a,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum", "AOS60:logDayNum"),
          order = c("minor", "AOS60", "AOA", "issues", "contig", "cincP", "polBinA", "nuk",
                    "nBellig", "keyAllies", "logBattles"),
          covariate.labels = c("Minor power", "Recent imbalance", "Overall imbalance",
                               "Issue salience", "Contiguity", "CINC ratio",
                               "Democratic initiator", "Nuclear", "Number of states",
                               "Major allies",
                               "Completed battles"),
          dep.var.labels.include=T,
          dep.var.labels = c("Third-party diplomatic intervention")
)